#include <csv.c>
#include <rnd/strconcat.c>
#include <muon1/particles.c>
#include <muon1/physics.c>
void report(const char *s) {printf("report: %s\n",s);}

void fixname(char *s)
{
	for (;*s;s++) if (!isalpha(*s) && !isdigit(*s)) *s='_';
}

typedef struct {char *name; double xoff,xp,zoff;} ElemParam;

LSet ElemParams_load(const char *filename)
{ // Get extra parameters from .lat (MAD format) file
	LSet ret=LSet(ElemParam);
	FILE *in=fopen(filename,"rb");
	char *p;
	for (char *s;s=fgetline(in);free(s))
		if ((stristr(s,"X_OFFSET") || stristr(s,"X_PITCH") || stristr(s,"Z_OFFSET")) && strchr(s,':'))
	{
		while (s[strlen(s)-1]==',') {char *t=fgetline(in); strconcat(&s,t); free(t);}
		ElemParam ep;
		ep.xoff=((p=stristr(s,"X_OFFSET"))?atof(p+strlen("X_OFFSET =")):0);
		ep.xp=((p=stristr(s,"X_PITCH"))?atof(p+strlen("X_PITCH =")):0);
		ep.zoff=((p=stristr(s,"Z_OFFSET"))?atof(p+strlen("Z_OFFSET =")):0);
		*strchr(s,':')=0; // Can't store this from if() because s might have changed
		streqi(&ep.name,s);
		fixname(ep.name); // Important that this is done here as well as in main
		LSet_add(&ret,&ep);
	}
	fclose(in);
	return ret;
}

#include <rnd/striseg.c>

double ElemParams_get(const LSet *es,const char *name,const char *param,const double d=0)
{
	ElemParam *e=(ElemParam *)es->a;
	for (int n=es->m-1;n>=0;n--) if (striseg(name,e[n].name))
	{
		if (!strcmpi(param,"X_OFFSET")) return e[n].xoff;
		if (!strcmpi(param,"X_PITCH")) return e[n].xp;
		if (!strcmpi(param,"Z_OFFSET")) return e[n].zoff;
		return d;
	}
	return d;
}

V3 getlayoutpos_bmad(csv l,const char *name)
{
	int n,m=csvm(l); char *c;
	for (n=1;n<m;n++)
	{
		fixname(l[n][0]);
		// layout_table (l) LA.CRMOD.GAT01 -> LA_CRMOD_GAT01
		// lat_table (name) LA.CRMOD.GAT01\1 -> LA_CRMOD_GAT01_1
		c=strrchr(name,'_');
		if (c && atoi(c+1) && !memcmp(name,l[n][0],c-name) || !strcmpi(name,l[n][0]))
			return V3_new(atof(l[n][2]),atof(l[n][3]),atof(l[n][4]));
	}
	return V3_new(-9e99,-9e99,-9e99);
}

#include <rnd/nonegi.c>

V3 getlayoutpos_muon1(csv l,const char *name)
{
	int n,m=csvm(l); char *t,*c;
	for (n=1;n<m;n++)
	{
		// Muon1 (l) LA_CRMOD_GAT01_1_E17
		// lat_table (name) LA.CRMOD.GAT01\1 -> LA_CRMOD_GAT01_1
		streqi(&t,l[n][1]);
		for (c=t+nonegi((int)strlen(t)-2);c>=t;c--) if (c[0]=='_' && c[1]=='E') {*c=0; break;}
		if (!strcmpi(t,name))
		{
			free(t);
			V3 x0=V3_new(atof(l[n][4]),atof(l[n][5]),atof(l[n][6])),
				x1=V3_new(atof(l[n][7]),atof(l[n][8]),atof(l[n][9]));
			double len=atof(l[n][2]),ang=atof(l[n][3]);
			if (len!=0 && ang!=0)
			{
				double R=len/ang,s=R-R*cos(0.5*ang);
				V3 p=V3_normalise(x1-x0)^V3_new(0,1,0);
				return 0.5*(x0+x1)+p*s;
			}
			else return 0.5*(x0+x1);
		}
		free(t);
	}
	return V3_new(-9e99,-9e99,-9e99);
}

#include <conio.c>

int main(void)
{
	ParticleSpecies_init();
	LSet extras=ElemParams_load("..\\cbeta_4pass_20160517.lat");
	csv f=csvload("..\\cbeta_4pass_20160517.lat_table",NULL),
		lbmad=csvload("..\\cbeta_4pass_20160517.layout_table"),
		lmuon=csvload("d:\\sbrooks\\muon1\\_Cbeta\\release20160517\\survey.csv");
	int n,m=csvm(f)-2;
	int cname=csvcol(f,"name"),ctype=csvcol(f,"key_name"),
		clen=csvcol(f,"L"),cdip=csvcol(f,"b_field"),cquad=csvcol(f,"b1_gradient"),
		cetot=csvcol(f,"e_tot"),
		clayx=csvcol(lbmad,"x");
	printf("%d elements\n",m);
	FILE *out=fopen("release20160517.txt","wt"),*outc=fopen("log.csv","wt");
	fputs("{Drift Length Angle}\nL90 0 90°\nR90 0 -90°\n",out);
	fputs("StartZ 15\nStartX 15.16\nStartAngle 0 -2.617994E-01\n",out);
	fputs("null\n\n",out);
	fprintf(outc,"n,Type,Name,Length (m),Dipole (T),Quad (T/m),BMAD x (m),y,z,\
Muon1 x (m),y,z,Etot (MeV),ke (MeV),p (MeV/c),kappa (rad/m),angle (deg),\
x_offset (m),z_offset (m),x_pitch (rad)\n");
	for (n=0;n<m;n++)
	{
		char *name=f[2+n][cname],*type=f[2+n][ctype]; fixname(name);
		double len=atof(f[2+n][clen]),b=atof(f[2+n][cdip]),quad=atof(f[2+n][cquad]);
		printf("\rConverting [%02d%%] %s",(int)(100.0*n/m),name); clreol();
		fprintf(outc,"%d,%s,%s,%.16lg,%.16lg,%.16lg",1+n,type,name,len,b,quad);
		V3 c=getlayoutpos_bmad(lbmad,name);
		if (c.x>-1e99) fprintf(outc,",%.16lg,%.16lg,%.16lg",c.x,c.y,c.z); else fputs(",,,",outc);
		c=getlayoutpos_muon1(lmuon,name);
		if (c.x>-1e99) fprintf(outc,",%.16lg,%.16lg,%.16lg",c.x,c.y,c.z); else fputs(",,,",outc);
		if (!strcmpi(type,"PATCH"))
		{
			double xoff,zoff,xp;
			fprintf(out,"%s_E%dX %.16lg\n",name,1+n,xoff=ElemParams_get(&extras,name,"x_offset"));
			fprintf(out,"%s_E%dZ %.16lg\n",name,1+n,zoff=ElemParams_get(&extras,name,"z_offset"));
			fprintf(out,"%s_E%dA 0 %.16lg\n",name,1+n,xp=ElemParams_get(&extras,name,"x_pitch"));
			fprintf(outc,",,,,,,%.16lg,%.16lg,%.16lg\n",xoff,zoff,xp);
		}
		else
		{
			double etot=atof(f[2+n][cetot])*eV,
				ke=etot-E0[Electron],
				p=p_from_ke(ke,E0[Electron]),
				kappa=charge[Electron]*-b/p,
				ang=kappa*len;
			fprintf(outc,",%.16lg,%.16lg,%.16lg,%.16lg,%.16lg\n",etot/MeV,ke/MeV,p/(MeV/C),kappa,ang radians);
			fprintf(out,"%s_E%d %.16lg %.16lg\n",name,1+n,len,ang);
		}
	}
	fclose(outc);
	fputc('\n',out);
	int patch,turn=0;
	fputs("Turn0: StartZ,L90,StartX,R90,StartAngle,\n",out);
	for (n=0;n<m;n++)
	{
		char *name=f[2+n][cname],*type=f[2+n][ctype];
		if (striseg(name,"LA_CRMOD_MAR_BEG")) {turn++; fprintf(out,"null;\nTurn%d: ",turn);}
		if (patch=!strcmpi(type,"PATCH"))
		{
			if (n%5!=0) fputc('\n',out);
			fputs("L90,",out); 
			fprintf(out,"%s_E%dX",name,1+n); 
			fputs(",R90,",out);
			fprintf(out,"%s_E%dZ",name,1+n);
			fprintf(out,",%s_E%dA",name,1+n); 
		}
		else fprintf(out,"%s_E%d",name,1+n);
		if (n+1<m) fputc(',',out); else fputc(';',out);
		if (n%5==4 || patch) fputc('\n',out);
	}
	fputs("\n\n",out);
	for (n=0;n<=turn;n++) fprintf(out,"Turn%d%c",n,(n<turn?',':';'));
	fclose(out);
	csvfree(f); csvfree(lbmad); csvfree(lmuon);
}